Advanced phase retrieval: maximum likelihood technique with sparse regularization of phase and amplitude

Sparse modeling is one of the efficient techniques for imaging that allows recovering lost information. In this paper, we present a novel iterative phase-retrieval algorithm using a sparse representation of the object amplitude and phase. The algorithm is derived in terms of a constrained maximum likelihood, where the wave field reconstruction is performed using a number of noisy intensity-only observations with a zero-mean additive Gaussian noise. The developed algorithm enables the optimal solution for the object wave field reconstruction. Our goal is an improvement of the reconstruction quality with respect to the conventional algorithms. Sparse regularization results in advanced reconstruction accuracy, and numerical simulations demonstrate significant enhancement of imaging. Introduction The conventional sensors detect only the intensity of the light, but the phase is systematically lost in measurements. Phase retrieval is a problem of the phase recovering using a number of intensity observations and some prior on the object wave field. The phase carries important information about an object shape what is necessary for a 3D object imaging and exploited in many areas such as microscopy, astronomy, etc. Moreover, phase-retrieval techniques are often simpler, cheaper and more robust comparing with interferometric ones. In 1982 Fienup introduced some, for now classical, iterative phase-retrieval algorithms [1]: error-reduction, gradient search and input-output methods. Many phase-retrieval methods are developed based on this pioneer work: the estimated magnitudes at the measurement planes are iteratively replaced by ones obtained from the intensity observations. * This work was supported by the Academy of Finland: project no. 213462, 2006-2011 (Finnish Programme for Centres of Excellence in Research) and project no. 138207, 2011-2014; and the postgraduate work of Artem Migukin is funded by Tampere Doctoral Programme in Information Science and Engineering (TISE). We are looking for the optimal wave field reconstruction from a number of intensity observations, and the reconstruction problem is formulated in terms of a variational constrained maximum likelihood (ML) approach. The spatial image resolution of the conventional phase-retrieval techniques is limited due to diffraction, what can be one of the main sources of artifacts and image degradations. In order to enhance the imaging quality and recover lost information, in this work we use the novel compressive sensing technique for the variational image reconstruction originated in [2]. The object wave field distribution is assumed to be sparse, and its amplitude and phase are separately analyzed and decomposed using very specific basis functions named BM3D-frames [3]. The proposed phase-retrieval algorithm is derived as a solution of the ML optimization problem using the BM3D-frame based sparse approximation of the object amplitude and phase distributions. Propagation model We consider a multi-plane wave field reconstruction scenario: a planar laser beam illuminates an object, and the result of the wave field propagation is detected on a sensor at different distances zr from the object at various measurement (sensor) planes. Here zr = z1+(r-1)· z, r=1,...K, where z1 is the distance from the object to the first measurement plane, z is the distance between the measurement planes, and K is a number of these planes. We assume that the wave field distributions at the object and sensor planes are pixel-wise invariant. In such a discrete-to-discrete model, the forward wave field propagation from the object to the r-th sensor plane is defined in the vector-matrix form as follows: 0 = , 1,... , r r r K u A u (1) where u0 and ur are n vectors, constructed by columns concatenating 2D discrete complexvalued wave field distributions (N×M matrices) at the object and sensor planes, respectively. Ar n×n is a discrete forward propagation operator, n=N M. We consider the paraxial approximation of the wave field propagation defined by the Rayleigh-Sommerfield integral. Depending of the used discretization of this integral, the operators Ar in (1) can be, for instance, the angular spectrum decomposition or the discrete diffraction transform in the matrix (M-DDT, [4]) or the Fourier transform domains (F-DDT, [5]). In our numerical experiments, we use F-DDT models enabling the exact pixel-to-pixel mapping u0 to ur. According to the used vector-matrix notation, the observation model with the additive zeromean Gaussian noise at the sensor planes, r[k] N(0, r2), takes the form: 2 =| | + , 1,... r r r r K o u (2) Let us assume that the object amplitude a0 n and phase 0 n can be separately approximated by small numbers of basic functions with coefficients a for the amplitude and coefficients for the phase. These basic functions are collected in the matrices a and for the amplitude and the phase, respectively. The amplitude and the phase are reconstructed from the noisy intensity data or. Decoupled augmented Lagrangain (D-AL) algorithm According to the maximum likelihood approach, the reconstruction of the object wave field is performed by minimization of the criterion: 2 2 2 a a 2 1 1 || | | || || || || || , 2 p p K r r l l r r L o u (3) subject to the following constraints: first of all the forward propagation models (1), and the constraints for the sparse modeling for the object amplitude and phase given as a0= a a, a= a a0, (4) 0= , = 0. (5) The quadratic fidelity term in (3) appears due to our assumption that the observation noise is Gaussian. The next two terms define the sparse regularization in the spectral domain, where the positive parameters a and define a balance between the fit of observations, the smoothness of the wave field reconstruction and the complexity of the solution. The first equations in (4) and (5) give the restrictions for the amplitude and the phase in the synthesis form and the second ones in the analysis form [2, 3]. Note that the priori unknown basic functions for the amplitude and the phase are selected from the synthesis a , and analysis a , matrices defining the redundant sets of the basic functions. The vectors a m are considered as spectra for parametric approximations of the amplitude and the phase. Conventionally, the sparsity of these approximations is characterized by a number of non-zero components of the spectra vectors (l0-norm, || ||0, of the vector ) or by the sum of the absolute values of the vector elements (l1-norm, || ||1 s| s|, of the vector ). In this work, the l1-norm is used based on the results stating that the solutions obtained for the l0or l1-norms are close to each other [6]. We use the augmented Lagrangian approach in order to reduce the constrained optimization for (3)-(5) to the unconstrained one. Furthermore, following the decoupling technique originated in [2, 3] instead of optimization of a single criterion we use the alternating optimization of two criteria L1 and L2: 2 2 2 2 1 2 r r 0 2 r r 0 0 0 2 2 1 1 1 1 2 1 = [ || | | || || || Re{ ( )}] || || 2 K H r r r r r r r L o u u A u u A u u v (6) 2 2 2 a a a a 0 2 r 0 2 a 1 1 || || || || || || || || 2 2 p p l l L (7) The criterion L1 in (6) is formed from the fidelity term and the forward propagation constraints similar to [7]. The variable v0 in the latter summand serves as a splitting variable separating optimization of L1 and L2 .This variable is an estimate of the object distribution u0. It is calculated as v0= a a exp(j ), where ‘ ’ stands for the element-wise multiplication of two vectors. , r, a and are positive parameters controlling “weight” of these penalty terms. ( ) stands for the Hermitian conjugate. A typical Lagrangian based optimization assumes: minimization of L2 with respect to a, ; minimization of L1 on u0, {ur}, and maximization of L1 with respect to the vectors of the Lagrange multipliers r n The splitting variable v0 separates minimizations on u0 and on a, . The solutions obtained for optimization of L1 and L2 results in the following iterative algorithm, similar to [2]: Initialize u0, { r} and calculate transform matrices a, , a, for t = 0 Repeat for t =1, 2...